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ABSTRACT 

In smooth-particle hydrodynamics (SPH), artificial viscosity is necessary for the correct treat- 
ment of shocks, but often generates unwanted dissipation away from shocks. We present a 
novel method of controlling the amount of artificial viscosity, which uses the total time deriva- 
tive of the velocity divergence as shock indicator and aims at completely eliminating viscosity 
away from shocks. We subject the new scheme to numerous tests and find that the method 
works at least as well as any previous technique in the strong-shock regime, but becomes vir- 
tually inviscid away from shocks, while still maintaining particle order. In particular sound 
waves or oscillations of gas spheres are hardly damped over many periods. 

Key words: hydrodynamics — methods: numerical — methods; A^-body simulations 



1 INTRODUCTION 

Smooth-particle hydrodynamics (SPH) is a Lagrangian method 
for mo dell ing fluid dyn amics, pioneered bv lGingold & MonaghanI 
( ll977l) and lLucvl ( ll977h . Instead of discretising the fluid quantities, 
such as density, velocity, and temperature, on a fixed grid as in Eu- 
lerian methods, the fluid is represented by a discrete set of moving 
particles acting as interpolation points. Due to its Lagrangian na- 
ture, SPH models regions of higher density with higher resolution 
with the ability to simulate large dynamic ranges. This makes it 
particularly useful in astrophysics, where it is used to model galaxy 
and star formation, stellar collisions, and accretion discs. 

The core of SPH is the kernel estimator: the fluid density is 
estimated from the masses m, and positions jc, of the particles vi^ 

p{xd = i:jmjW{\x,-Xj\,h,\ (1) 

where W is the kernel function and /j, the SPH smoothing lengtl[| 
for the ith particle. Similar estimates F{x) for the value of any field 
Fix) can be obtained from its discretised values F,. By applying 
these estimators to the fluid equations governing mass, momentum 
and energy, discrete equations for the SPH particle positions x, and 
other properties (such as internal energy k,) can be obtained. To- 
gether with an appropriate time integration method, these constitute 
a concrete SPH scheme. 

Unfortunately, this process is not unique and since its incep- 
tion the SPH method has undergone many refinements such as indi- 
vidual particle smoothing lengths and viscosities, as well as many 
alternative derivations of the SPH equations, leading to a plethora 
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' We use the symbol " to denote a local estimate - in many SPH-related 
publications the distinction between actual and estimated quantities is not 
clearly made, confusing the discussion. 

^ In this study we use the convention that the kernel has finite support of 
one smoothing length radius, i.e. W = for \xi - Xj\ > h. 
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Figure 1. A ID sinusoidal sound wave with velocity amplitude lO^^c and 
y = 1.4 propagated for 50 periods with SPH without artificial viscosity 
using 100 particles and with a grid code (Ramses, Tevssier 2002) using 128 
cells (only every fifth particle or grid cell is plotted). Both methods preserve 
the wave amplitude and period, demonstrating their dissipation-less nature. 



of SPH methods. While formally these various schemes differ only 
in their error terms, their conservation and stability properties can 
be quite different. This has lead to the unfortunate situation that the 
shortcomings of a few such implementations are often blamed on 
the general SPH concept p er se. 

However, iSpringel & Hemguistl ( |2002|) have pointed out that 
SPH equations derived from a variational principle are not only 
unique, but also conservative. Such SPH equations are most simply 
obtained as the Euler-Lagrange equations derived from an SPH La- 
grangian £. representing the Lagrangian of the fluid system. Once 
£. is chosen, the SPH equations follow uniquely (see Appendix I A2 1 
for a typical example). Complementing these with a symplectic in- 
tegrator, such as the standard leap-frog, results in a SPH scheme 
which by construction conserves the total mass, momentum, angu- 
lar momentum, energy, and entropy. 

The conservation of entropy means that SPH is dissipation- 
less, as demonstrated in Fig. [T] In real fluids, however, entropy in- 
creases in shocks, where particle collisions randomise their veloci- 
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ties generating hieat and entropy. Tiiis basic collisional meclianism 
is inlierent to all fluids (except for dust and collisionless plasma, 
which therefore may not be considered fluids) and prevents the flow 
from becoming multi-valued. In SPH artificial viscosity is needed 
to dissipate local velocity differences and convert them into heat, 
which generates entropy and prevents inter-penetration of SPH par- 
ticles and thus a multi-valued flow. 

Since the artificial viscosity required for this goal is usually 
much stronger than the actual physical viscosity, it also causes un- 
physical dissipation away from shocks. While it may be possible 
for certain simulations to select the magnitude of the viscosity to 
minimise such undesired dissipation, in general the adverse effect 
of artificial viscosity is unknown prior to any simulation and, pos- 
sibly, even afterwards. For example, when simulating the effect of a 
perturbing massive body on a pulsating star, it may be very difficult 
to distinguish this effect from that induced by artificial viscosity. 
Another example is the case of a differentially rotating disc, where 
artificial viscosity causes spurious angular momentum transport. 

Since viscosity is a dissipative process, the corresponding 
SPH equations cannot be derived from a variational principle, and 
we are back to ad-hoc methods for deriving them. Most SPH simu- 
lations to date still use a rather simple artificial viscosity, which ef- 
fectively amounts to modelling a viscous fluid and quickly damps 
away any oscillations, such as sound waves or stellar pulsations, 
and impedes shear flows. While suggestions have been made to re- 
duce such unwanted dissipation, our goal here is to eliminate it. To 
this end we introduce a novel method of controlling the amount of 
artificial viscosity, such that away from shocks the modelled flow 
is virtually inviscid. 

Section |2] describes SPH artificial viscosity and previous ef- 
forts to reduce its adverse effects, while our new method is outlined 
in Section[3] The ability of the new scheme to reduce artificial vis- 
cosity but also to capture shocks is demonstrated in Sections |4]and 
[S] respectively. Finally, Section[6]concludes our study. 



2 REDUCING UNWANTED ARTIFICIAL VISCOSITY 
2.1 Standard SPH artificial viscosity 

The traditional form of artificial viscosity (e.g. lMonaghaijri992t) 
adds the following terms to the momentum and energy equations, 
allowing the conversion of kinetic energy into heat. 



('>;)av = -H^m^Hi, V, lVy 

(".■)aV = 7 Yjj mj Hy Vij ■ V, VK,y 

with the average kernel 

W^j= {{W{\x,j\,h,) + W(\x.j\,hj)). 



(2a) 
(2b) 

(3) 



Here, Xy = jc, - Xj and Vy = w,- -Vj, while /i,- is the individual adap- 
tive smoothing length of each SPH particle (for details on how hj is 
adapted see AppendixlAU. T he artificial viscosity term is given by 
dGingold & Monaghaij|l983h 



with 



for Vij ■ Xjj < 



otherwise 



hjj Vjj ■ Xjj 
xl + e2 



(4) 



(5) 



(Jiij = [hi + hj]/2 and likewise for the average sound speed c,y and 
estimated density py). Since Tly = for receding particle pairs, ar- 
tificial viscosity does not affect expanding flows. This functional 
form of SPH artificial viscosity may seem rather ad-hoc, but it 
is reasonably well motiva ted and emerged as the mos t useful one 
amongst several methods dOingold & Monaghanlll983i) . Moreover, 
it is eq uivalent to the f orm of dissipation implicit in Riemann 
solvers (Monaghanl l 19971) . 

By expanding density and velocity in a Taylor series around 
Xj, it is straightforward to show that these terms correspond to both 
a shear and a bulk viscosity. More quantitatively, if one assumes 
that, other than in equation l|4j, artificial viscosity acts between ap- 
proaching and receding neighbours and that /? = 0, t he correspond- 
ing s hear and bulk viscosity coefficients are (e.g. iMeglicki et al.l 
1 19931) 77 = \aKhcp and (, = ^r}, respectively, where the factor k is of 
order unity and depends on the functional form of the kernel. This 
implies that artificial viscosity decreases with increasing resolution 
(smaller h). Thus, a straightforward though expensive way to re- 
duce unwanted dissipation is to increase the resolution. In fact, one 
motivation for reducing artificial viscosity is to avoid this purely 
numerical necessity for high resolution. 

Most SPH applications to date us e the above trea tment with 
a = I. The widely used code gadget-2 ( ISpringelll2005|l employs a 
fixed a chosen at the start of the simulation (though llDolag et ^ 
1 20051 have implemented into gadget-2 the improved method de- 
scribed in 32.3l below). Clearly, in complex situations, where strong 
and weak shocks are present as well as converging flows, any 
choice for a is unsatisfactory, leading to bad treatment of strong 
shocks, over-damping of converging flows, or both. 



2.2 Balsara's method 

The purpose of artificial viscosity is to allow for entropy gener- 
ation across shocks and to stop particle interpenetration. To this 
end, only bulk viscosity is required, but the inherent shear viscosity 
is unnecessary. What is worse, this shear viscosity may seriously 
compromise simulations of shear flows, such as in a differentially 
rotating gas disc. In an eff ort to reduce the resulting artificial angu- 
lar momentum dissipation. ISalsaral ( Il995h proposed to multiply H,j 
with a reduction factor /, = (/ -I- fj)l2 with 



f. 



|V-y,| + |Vxw,| 



(6) 



(with velocity divergence and curl estimated using the SPH ker- 
nel estimator). This term diminishes the effect of artificial viscosity 
whenever the vorticity dominates the convergence. However, this 
method only reduces (but does not eliminate) unwanted dissipation 
in the presence of a rotating shear flow. 



2.3 Tiie method of iMorris & Monaghanl 

Standard SPH artificial viscosity acts whenever the flow of the fluid 
converges, even if only weakly. For example, when a pulsating star 
contracts artificial viscosity damps its pulsation. Exactly the same 
happens to ordinary sound waves: standard SPH viscosity damps 
them, as demonstrated in Fig. |2l the faster the shorter the wave 

length (because these ar e more poorly resolved ). 

With this in mind, iMorris & Monagtianl j 19971) proposed to 
adapt the strength of artificial viscosity to the local convergence 
of the flow. To this end, they introduced the concept of individual 
adaptive viscosities o-, for each particle, replaced a in equation Q 
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Figure 2. As Fig.Q] but for SPH with standard (a = 1) or Moms & Mon- 
aghan (1997) artificial viscosity, as well as our new method (only every fifth 
particle is plotted). Also shown are the undamped wave (solid) and lower- 
amplitude sinusoidals (dashed). Only with our method the wave propagates 
undamped, very much like SPH without any viscosity, as in Fig.[T] 



by Q-y = (a, + aj)/2, and set /3 oc o-y . The individual viscosities are 
adapted according to the differential equation 



ai = (ffmin - + Si 

with the velocity-based source term 
S, = maxj - V-w„ 0). 
and the decay tim^l 

T, = h,/{2ec,). 



(V) 



(8) 



(9) 



Here, = 0. 1 constitutes a lower limit for the artificial viscos- 
ity such that or, = (jr„,jn for non-convergent flows. For a convergent 
flow, on the other hand, a; grows above that value, guaranteeing the 
proper treatment of shocks. In the post-shock region, the flow is no 
longer convergent and ff, decays back to on the time scale t, 
(typically f = OA - 0.2). This method reduces the artificial viscos- 
ity away from shocks by an order of magnitude compared to stan- 
dard SPH and gives equally accurate post and pre-shock solutions 

jMorris & Mon aghan 199J). 

More recently, Rosswog. Davies. Thielemann & Piraij ilOOd) 
proposed to alter the adaption equation lO tc0 



min 

)/Ti + (a 

max 

1.5, while lPricd ^20041) advocated < 



(10) 



= 2. The effect 
max and second to 



tt, = {a, 
with a, 

of this alteration is first to prevent a, to exceed a 
increase a, for small a,, which ensures a faster viscosity growth, 
resulting in somewhat better treatment of shocks ( Price 2004). This 
method may also be combined with the Balsara switch by applying 
the reduction factor (|p either to H,, (Rosswog e t al., 2000.) or to S,- 
dMorris & MonaghanlIl997l : IWetzstei n et al. 20(^7 



The scheme of equations {SJ, (O and dlOl l with aj„i„ = 0.1, 
Q^max = 2 and f = 0.1 is the current state of the art for SPH and 
is implemented in the codes phantom (by Daniel Price) and vine 
dWetzstein et alj2009i) . In sections|4]and|5l we will frequently com- 
pare our novel scheme (to be described b elow) with t his method and 
refer to it as the 'M&M method' or the ' Pried ( |2004|) version of the 
M&M method' as opposed to the 'original M&M method', which 
uses equation ([7} instead of dlOt . 

^ The factor 2 in the denominator of equation j9) accounts for the dif- 
ference in the definition of the smoothing length h between us and 



iMorris &Monaghai]h997h . 
^ This is equivalent to k eeping j?) but m ultiplying the source term (S) by 
(fmax ~ <^), which is what lRosswog et al] actually did. 



2.4 Critique of the M&M method 

The M&M method certainly constitutes a large improvement over 
standard SPH, but low- viscosity flows, typical for many astrophys- 
ical fluids, are still inadequately modelled. After studying this and 
related methods in detail, we identify the following problems. 

First, any a„^m > results in unwanted dissipation, for example 
of sound waves (see Fig.O or stellar pulsations (see yet the 

M&M method requires o-min « 0. 1 . This necessity ha s been estab- 
lished by numerous tests (most notablvof |Pricel2004l) and is under- 
stood to originate from the requi rement to 'maintain order am ongst 
the particles away from shocks ' (^ orris &Monaghaj|l997l) . 

Second, there is a delay between the peak in the viscosity a 
and the shock front (see Fig. [3}: the particle viscosities are still 
rising when the shock arrives. One reason for this lag is that inte- 
grating the differential equation dlOl l increases a, too slowly: the 
asymptotic value 



1 SiTi 



1 + SiTi 



(11) 



is hardly ever reached before the shock arrives (and 5; decreases). 

Third, the source term ([8j does not distinguish between pre- 
and post-shock regions: for a symmetrically smoothed shock it 
peaks at the exact shock position (in p ractice the peak occurs one 
particle separation in front of the shock. Morris & Monagh anI 19971 . 
see also Fig.[3j. However, immediately behind the shock (or more 
precisely the minimum of V-w), the (smoothed) flow is still con- 
verging and hence a continues to increase without need. A further 
problem is the inability of the source term ^ to distinguish be- 
tween velocity discontinuities and convergent flows. 

Finally, in strong shear flows the estimation of the velocity di- 
vergence V-w, needed in ([8}, often suffers from substantial errors 
(see Appendix IB II for the reason), driving artificial viscosity with- 
out need. This especially compromises simulations of differentially 
rotating discs even when using the Balsara switch. 



3 A NOVEL ARTIFICIAL VISCOSITY SCHEME 

Our aim is a method which overcomes all the issues identified in 
32.4l above and in particular gives a, — > away from shocks. To this 
end, we introduce a new shock indicator in 33.11 a novel technique 
for adapting a, in 33.21 and a method to suppress false compression 
detections due to the presence of strong shear in 33.31 



3.1 A novel shock indicator 

We need a shock indicator which not only distinguishes shocks 
from convergent flows, but, unlike V-y, also discriminates between 
pre- and post-shock regions. This requires (at least) a second-order 
derivative of the flow velocity and we found the total time deriva- 
tive of the velocity divergence, V v = d(V-i')/df. to be most useful. 
As is evident from differentiating the continuity equation. 



(12) 



V v < indicates an non-linear density increase and a steepen- 
ing of the flow convergence, as is typical for any pre-shock region. 
Conversely, in the post-shock region V v > 0. This suggests to 
consider only negative values and, in analogy with equation {Sj, we 
define the new shock indicator 



Ai = max j - V vi, 0|. 



(13) 
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Here, ^, is a limiter, detailed in 93.3 1 below, aimed at suppressing 
false detections of compressive flows in multi-dimensional flows. 



3.2 Adapting individual viscosities 

Instead of increasing a, by integrating a differential equation, we 
set a, directly to an appropriate local value o-ioc,; whenever this ex- 
ceeds the current value for or,. After extensive experimenting, we 
settled on the following simple functional form 



'u'- . + h^A, 

sig.f ( ' 



with the signal velocit}!! 

Wsig.j = rnax (c„ - minjO, Vy ■ Xy} 



(14) 



(15) 



At the moment of passing through a shock (more precisely through 
a maximum of the flow convergence), A and hence aioc return to 
zero and whenever a, > a\acj we let a, decay according to 



a/ = (ffloc./ - «/)/•!",■, 



(16) 



We use fsig, rather than c in the decay time t, for internal con- 
sistency (this is of little practical relevance as v^\^ x c in the post- 
shock region). We use C = 0.05, such that the viscosity decays twice 
as slowly as in previous methods, avoiding some occasional minor 
post-shock ringing not present in methods with a^i„ > 0. However, 
the traditional f = 0.1 also gives satisfactory results for most of our 
test problems. 



3.3 Avoiding false compression detections 

As explained in detail in Appendix IBll in multi-dimensional flows 
strong shear induces false detections of V v with the standard SPH 
estimator even in the absence of particle disorder (noise). As shown 
in Appendix |B2l these errors can be reduced by first estimating the 
velocity gradient matrix V = Vv and then obtaining V-w as its trace 
(we employ a similar method to estimate V-v, see Appendices lB3b . 

Unfortunately, even with this improved method false detec- 
tions for V v (and V-w) remain, for example in the situation of a 
differentially rotating disc. These still induce artificial viscosity, 
which may be significant in particular if cjh is small compared 
to the shear. The limiter ^, in equation l |13t is aimed at suppress- 
ing such false detections by ^, — > whenever the shear is much 
stronger than the convergence and no shock is present. 

Having obtained the velocity gradient matrix V, the shear is 
easily obtained as its traceless symmetric part S = (V + V')/2 - 
V"' (V-w)l (with V the number of spatial dimensions), while the pres- 



ence of a shock is indicated by 



-1 « ij, = 1 y sign(V-i;j) mj WQx, - Xj\, h,}. 



(17) 



since near a shock V i/ < for all particles. After some experiment- 
ing, we found the following functional form for the limiter suitable 



\2(\-R,fV-v,\^ + lr{SrS\) 



(18) 



^ Vaiious definitions for the signal velocity can be found in the SPH liter- 
ature. Ours reflects the maximum velocity with which information can be 
transported between particles, but avoids i/jig , ^ 0. 



This functional form is similar to the Balsara limiter ^ in that it 
compares the flow convergence to a measure of the traceless part of 
the velocity gradient (the shear or the vorticity). 

Alternatively, if one can be sure that no strong shear flows oc- 
cur during the simulation, one may use the standard SPH estimator 
for V-y and estimate V-u from its change over the last time step. 
However, the limiter is still desirable and one may use |V x in- 
stead of tr(S-S') in equation l |18t . We do not use this simplified 
version in the tests presented below, but our experiments indicate 
that such a method would pass all our tests except those of 94. 3 l and 
3] both involving strong shear. 



3.4 Behaviour in typical situations 

Before considering 2D and 3D test problems, we now assess the 
behaviour of our novel scheme, as well as that of the M&M method, 
in simple yet typical situations. 

First, consider a sound wave of velocity amplitude <c c 
and wave number k <c as example of a well-resolved weakly 
convergent flow. In this case, A ^ k'cvs and 5 - kv,, at their respec- 
tive maxima. Since u^\^ ^ c » we have a]„^ ^ ay„.^jfk^(ujc), 
while for the M&M method the asymptotic value - ffmin + 
Q'maxhk{vJc)/2{. In the limit kh — » of a well-resolved wave, 
Q'loc faster than — > ffmin, such that even with a^i„ = 
the M&M method would be more viscous than our new scheme. 
Fig.[2]shows ID sound-wave SPH runs, demonstrating that our new 
scheme behave s quasi-inviscid in thi s situat ion. 

Following lMorris & MonaghanI ( Il997l) , we may also consider 
a simple homologous flow v = -ax with a > 0, which approxi- 
mates certain astrophysical problems involving collapse and does 
not require artificial dissipation. For this situation S =3a but A = 
(a direct consequence of the ability of V-y to distinguish shocks 
from convergent flows), such that our new scheme remains invis- 
cid, while the M&M method does not even for ay„i„ = 0. 

Next, consider a strong shock with velocity discontinuity 
Sv » c. Assuming that it is smoothed over one kernel width, we 
find maximum amplitudes of S - 5v/h and A ^ {Svjhf' (the exact 
values depend on the shock conditions and the functional form of 
the smoothing kernel). Since Vs\^ ^ hV-v ~ 6v, our new scheme 
gives Q'loc ~ Qmax, whilc the asymptotic value i ll It for the M&M 
method also approaches o-max- 

While 3D simulations of strong shocks are presented in 95.21 
Fig. [3] presents weak ram-shock simulations with Su = 0.1c (top) 
and Su = c (bottom) for our new scheme, the M&M method, and 
standard SPH. In both regimes the peak in, respectively, ajoc and 
a, is one particle farther in front of the shock with our new method 
than with the M&M method, which reflects the superiority of A 
over S to detect an incoming shock. This, combined with setting 
the viscosity parameter directly to the required value, results in the 
peak in a to occur two particle separations before the shock for our 
new method, while for the M&M method it peaks a similar length 
behind the shock. 

With our new method, the viscous deceleration (bottom pan- 
els in Fig.[3]l sets in about three particle separations before the weak 
and the strong shock, yielding good shock capturing properties in 
both cases. The M&M method, on the other hand, decelerates the 
flow much earlier for a weak shock than for a strong shock and re- 
sults in significant over-damping of weak shocks (which also per- 
tains to density and internal energy - not shown in Fig. [3]l, while 
our method smoothes both shocks over four particle separations 
(top panels in Fig. [3), the optimal SPH resolution in ID. Note that 
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Figure 3. A ID ram shock with Sv = Q.l (top) and 6v = 1 (bottom) in 
ideal gas with y = 1.4 simulated using standard SPH, the M&M and our 
new method. We compare the velocity, viscosity parameter, its asymptotic 
value and the viscous deceleration. Initially, the velocity is discontinuous 
with V = -Sv sign(ji:), resulting in two shocks of Su propagating in either 
direction from the origin; the shock plotted propagates from right to left. 

standard SPH is hopeless: it over-smoothes the strong shock and is 
completely incapable of dealing with the weak shock. 

3.5 Maintaining particle order 

The main point of our method is the absence of artificial viscosity 
away from shocks. Hence, if o-n^in > was indeed required to main- 
tain particle order, as previously argued in context of the M&M 
method, our method should fail in this regard. Noise in SPH can 
emerge from shocks or carelessly generated initial conditions. 

Let us first consider the time evolution of noisy initial condi- 
tions, generated by adding random displacements to particle posi- 
tions representing noise-free hydrostatic equilibrium (the vertices 
of a face-centred-cubic grid, i.e. densest-sphere packing). We con- 
sider two cases with the displacements in each dimension drawn 
from a normal distribution with rms amplitude equal to the nearest- 
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Figure 4. Time evolution of </mi„, defined in equation 09) , for SPH simula- 
tions started from noisy initial conditions (see text). All SPH schemes with 
ailificial viscosity suppress the noise equally well. 
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Figure 5. The rms amplitudes of density and velocity fluctuations for 3D 
simulations of the Sod (1978) shock tube test (see also Fig. II 11 . Initial con- 
ditions were prepared using a glass. The shock propagates to the right and 
is indicated by the dotted line; the velocity jump at the shock is 0.63. 



neighbour distance or a tenth of it, respectively. The time evolution 
of such noisy initial conditions can be distinguished by monitoring 

q„,i„ = mm{\Xij\/hi]. (19) 

'J 

There are three possible scenarios. Either the particles settle back 
close to the original grid (qmin approaches its grid value ^grid), form 
a glass (^niin approaches a finite value < ^gnd), or form dense 
clumps ('clumping instability', q^;„ ~ 0). Fig. |4]plots the evolu- 
tion of ^niin for Nf, = 40 SPH neighbours (see also Appendix I A It 
when ^grid ~ 0.529. Clumping only occurs when a = 0, while for 
any viscous scheme tested the particles settle back onto the grid or 
form a glass with roughly similar time evolutions. 

Post-shock noise occurs because the shock-induced compres- 
sion disrupts the original particle order, but other than in the above 
test the viscosity is already switched on. In Fig. |5] we plot the am- 
plitudes of the velo city and density noise in 3D simulations of the 
standard fSo j d 1 97§) shock tube test (see also 95.lt . The three meth- 
ods have similar levels of density noise, but standard SPH is less 
noisy in the velocities, which is not surprising given its stronger 
viscosity. However, between the two viscosity suppressing schemes 
there is little difference, even though Q-„,i„ = for our method. Sim- 
ilar results obtain for other shock tests and we conclude that our 
method is no worse than M&M's for maintaining particle order. 
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Figure 6. Steepening of a ID sound wave: velocity and viscosity param- 
eter vs. position for standard SPH, the M&M method, our new scheme, 
and Godunov p ailicle hydrodynamics of first and second order (GPH, 
ICha & Whitworth 2003), each using 100 pailicles per wavelength. The soUd 
curve in the top panel is the solution obtained with a high-resolution grid 
code. 



4 VISCOSITY SUPPRESSION TESTS 

We now present some tests of low-Mach-number flows, where pre- 
vious methods give too much unwanted dissipation. 



4.1 Sound-wave steepening 

The steepening of sound waves is a simple example demonstrat- 
ing the importance of distinguishing between converging flows and 
shocks. As the wave propagates, adiabatic density and pressure os- 
cillations result in variations of the sound speed, such that the den- 
sity peak of the wave travels faster than the trough, eventually try- 
ing to overtake it and forming a shock. 

In our test, a ID sound wave with a velocity amplitude 10% of 
the sound speed is used (ideal gas with y = 1.4). Fig.|6]compares 
the velocity field at the moment of wave steepening for various SPH 
schemes, each using 100 particles, with a high-resolution grid sim- 
ulation. The new method resolves the shock better than the Mc&M 
scheme, let alone standard SPH. 

In Fig. [6] we als o show results from GPH (Godunov- type par- 
ticle hydrodynamics, ICha & Whitworthll2003h , which differs from 
SPH by using the pressure P', found by solving the Riemann prob- 
lem between particle neighbours, in the momentum and energy 
equations and avoids the need for explicit artificial viscosity. This 
subs titution does not affect the energy or momentum conservation 
jCha 2002), and indeed we find that both are well conserved. While 
the first-order GPH scheme is comparable to standard SPH and also 
to an Eulerian Godunov grid code using the same Riemann solver 
without interpolation (not shown), the second-order GPH scheme 
resolves the discontinuity almost as well as our novel method. 



4.2 ID converging flow test 

Similar to sound-wave steepening, this test requires good treatment 
of convergent flows and weak shocks. The initial conditions are 
uniform pressure and density and a continuous flow velocity 



4{l+x)v„ -1.00 < X < -0.75, 

Va -0.75 < X < -0.25, 

-4xv„ -0.25 < X < 0.25, (20) 

-y„ 0.25 < X < 0.75, 

4(l-x)v„ 0.75 <x< 1.00. 
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Figure 7. A ID converging flow test with initially constant density and 
pressure and velocities given by equation )20t using an adiabatic equation 
of state with y = 1.4. Top: run for v^, = 1 at ( = 0.3; bottom: run for 
= 2 at ? = 0. 1 . The solid lines are the result of a high-resolution Eulerian 
grid-code simulation. 



As there is no analytical solution, we compare the results to a high- 
resolution grid-code simulation. We run tests for v„ = I and v„ = 2 
as shown in the top and bottom panels of Fig. [7] 

While the M&M switch certainly improves upon standard 
SPH, it still over-smoothes the velocity profile as the viscosity is 
increased before a shock has formed. This is particularly evident in 
the velocity profile of the v^, = 2 case (bottom) near x = 0. The 
new switch keeps the viscosity low, in the u„ = 2 case an order 
of magnitude lower than the M&M method. In fact, the agreement 
between our method and the high-resolution grid code is as good 
as one can possibly expect at the given resolution, in particular the 
velocity plateau and density amplitude around x = in the = 2 
case (bottom) are correctly modelled. 
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Figure 8. Keplerian ring test: particle positions at various times for standard SPH with Balsara switch, the M&M method with and without Balsara switch, 
and our new method without and with the viscosity limiter f of equation (18). Only for this last method the ring remains stable against a viscosity-induced 
instability. (Ring-like features at r < 2 are artifacts caused by the dynamical time close to the centre being short compared to the time step). 



4.3 2D Keplerian-ring test 

In this test, a gaseous ring orbits a central point mass, neglect- 
ing the self-gravity of the gas. Initially, the ring is in equilib- 
rium: pressure forces, attraction by the point mass, and centrifu- 
gal forces balance each other. The Keplerian differential rotation 
implies that t he flow is shearing and any viscosity causes the 
ring to spread (Lvnden-Bell & Pringle 1974). This is indeed what 
[Maddison, Murray & M onaghan. (.1996.) found in SPH simulations 
without pressure forces. 

iMaddison et al.l also found an instability to develop 
from the inner edge, which quickly breaks up the ring. 
They argue convincingly that this is the viscous instability 
dLvubarskii. Postnov & Prokhorovl[T994h . which causes eccentric 
orbits at the inner edge of the ring to become more eccentric due to 
the viscous deceleration peaking at apo-centre. 

Ilmaeda & Inutsukal ( |2002|) performed SPH simulations of the 
same problem but including pressure forces. They find a similar 
break-up of the ring after only few rotations and blame it on an 
inadequacy of the SPH scheme itself. We strongly suspect that 
Ilmaeda & Inutsukal encountered a form of the clumping instabil- 
ity, which appears to be particularly strong in 2D simulations of 
strong shear flows (though it may have been a d ynam ical instabil- 
ity inherent to gaseous Keplerian rings, e.g. [Papaloizou & Pringld 
ll984ll985l : [Goldreich & Naravanlll985h . This numerical instabil- 
ity grows on a local hydrodynamical time and may therefore be sup- 
pressed by choosing the sound speed c much lower than the rota- 
tion speed u^. Indeed, Price (2004) and Monaghan (2006), who re- 
peated these and similar experiments with a very low sound speed, 
found no such numerical instabilities. A detailed investigation of 
these issues is clearly beyond the scope of our study and we merely 
compare our new scheme to previous methods for pressure forces 



with c <c when the viscous instability should strike after few 
rotations depending on the strength of the artificial viscosity. 

In our test, GM = 1000 for the central point mass, 
while the gas ring has Gaussian surface density centred on 
r = 10 with width (standard deviation) 2.5 represented by 
= 9745 particles initially placed accord ing to the method of 
'Cartwright. Stamatel los & Whitwortlj (12009'). This implies an or- 
bital period of 2;r and velocity of = 10 at the ring centre. We 
choose a sound speed of c = 0.01 <c to ensure that any dynami- 
cal instabilities of inviscid rings become important only after many 
periods. 

Figure [8] shows the particle distributions at various times for 
different SPH schemes. Only with our new method, the rings stay 
in their initial equilibrium configuration over at least five periods, 
while for the other methods, the inner parts of the ring soon become 
disordered leading to a catastrophic break-up after a few periods. It 
is noteworthy that this failure occurs despite the Balsara switch, 
which was designed specifically for applications like this. 

Note that without the viscosity limiter £ of equation l |18l l, our 
novel method fails, precisely because of shear causing false detec- 
tions of flow compression (as highlighted in 93.3l and Appendix iBt. 

We also run similar tests with the central point mass replaced 
by a mass distribution (Plummer sphere or Kuzmin disc) with grav- 
itational potential O = -GM/ y/r^ + with s = 3, such that the 
rotation curve of the disc also contains a rising part, similar to the 
situation in galactic discs. The outcome of these simulations (not 
shown) is essentially identical to that for the pure Keplerian rings: 
only our new method with viscosity limiter does not fall prey to the 
viscous instability. 
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Figure 9. Left: virial ratio plotted versus time for SPH models of a radially oscillating polytrope which initially was in its fundamental radial eigenmode with 
relative radial amplitude of 0.01 and period 3.89. The solid curves are for a SPH model without any artificial viscosity. Right: the viscosity parameter a at 
t = 91 (maximum contraction) and t = 99 (maximum expansion) for every 100th particle. The new method keeps viscosity lower at the edge of the polytrope. 




Figure 10. Same as Fig.|9] except that the sphere is in circular' orbit around a point mass of 100 times its mass and with orbital radius 20 times its radius (the 
kinetic and potential energies are connected for the contributions from the orbit). The viscosity parameter for every 100th particle is plotted alt = 100 (riglit). 



4.4 An oscillating polytropic sphere 

The pulsations of a polytropic s phere are a good test for the adverse 
effects of artificial dissipation iSteinmetz & Mtilieilll993l) . We set 
up a polytropic sphere o f 10^ partic les and induce oscillations in its 
fundamental mode (e.g. ICoxIl 19801) with relative amplitude of 0.01 
in radius and a period of P = 3.8. 

In the absence of viscosity we expect the radial oscillations 
to continue with the initial amplitude and period over many os- 
cillations. However, as with any numerical method some small 
amount of numerical dissipation may appear. Nonetheless, such ef- 
fects should be small compared to the dissipation caused by arti- 
ficial viscosity. Since the size of the radial perturbations increases 
with radius, we expect the oscillations to be small at the centre of 
the polytrope and therefore our new method to keep the viscosity 
low there. However, at the edge the size of the oscillations is more 
significant, and we may see an increase in viscosity at this point. 



In order to track the oscillations, we monitor in Fig.|9]the time 
evolution of the virial ratio -2{T + U)IW where T, U, and W, are 
the kinetic, internal, and the gravitational energies, respectively. At 
maximum contraction the virial ratio is at its peak and at maxi- 
mum expansion the virial ratio is lowest. With no artificial viscosity 
(solid curves in Fig. [9} the wave remains at constant amplitude bar- 
ring a slight initial drop. The period averaged over 25 oscillations 
is P = 3.89, only slightly larger than the expected value. The rea- 
son for this discrepancy is most likely the unavoidable deviation of 
the (finite-resolution) SPH model from a perfect polytropic sphere. 
This deviation also means that our SPH model is not exhibiting a 
pure eigenmode, but in addition contains some higher-order modes 
at low amplitudes, resulting in some beating between them. 

The M&M method results in a slow but continuous decay of 
the oscillations, though the period is hardly affected. This damping 
can be blamed largely on the finite ttmin (standard SPH damps the 
oscillation ten times faster). Conversely, our new method, hardly 
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Figure 11. Comparison of our new scheme and the M&M method for the standard Sod (1978) shock tube test with the analytic solution (solid). 



damps the oscillations at all, because a is kept very small (except 
for the outermost layers where a is still below the M&M values). 

We also run simulations where the oscillating polytropic 
sphere is on a circular orbit 20 times the radius of the sphere 
around a point mass 100 times that of the sphere (corresponding 
to a period of 56 time units). With this choice, the tidal radius is ap- 
proximately four times the radius of the gas sphere, implying that 
tides are strong but not catastrophic. Since the orbital accelerations 
are much larger than those due to the polytropic oscillations, this 
is a tough test for any numerical scheme. In particular, Eulerian 
methods should have severe problems (this does exclude using co- 
rotating coordinates, which do not allow for tidal evolution of the 
orbit and are unavailable for eccentric orbits). 

The time evolution of the virial ratio and the viscosity param- 
eter a are shown in Fig. [TO] for the same viscosity schemes as for 
the isolated case in Fig.|9] First note that the undamped simulations 
(solid curves) behave differently from the isolated case, exhibiting 
variations and a slight decay, both of which are most likely caused 
by the tidal field. As to be expected for any Lagrangian scheme, 
both SPH methods perform very similar to the isolated case, be- 
cause neither V-u nor V-y are afi'ected by the orbital acceleration. 
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Figure 12. Same as Fig. [5] but for M = 50. We distinguish between the 
original M&M method (using eq. |7) and the Price (2004) version (using 
eq. llOl with ttmax = 2), which has been denoted 'M&M' in all figures so far. 



5 SHOCK CAPTURING TESTS 

In this section, we subject our method to situations where artificial 
viscosity is required, mainly high-Mach number shocks, and our 
aim is to demonstrate that it performs at least as well as the M&M 
method. 

5.1 Sod shock tube test 

Theiol ( Il978h shock tube test is a standard test for any shock 
capturing method and consists of an initial discontinuity in pressure 
and density leading to the production of a rarefaction wave, contact 
discontinuity and shock wave, which forms from the steepening of 
a subsonic wave. The whole system is subsonic with a maximum 
Mach number of M ~ 0.63 in the pre-shock region. We perform the 
test in 3D at a resolution of 200 particle layers in the high-density 
region. 



The density, energy, velocity, and viscosity for standard SPH 
as well as the M&M and our method are shown in Fig. [TT] As for 
the ID ram test (see FigO, our new method switches on viscos- 
ity already in the pre-shock region peaking about one smoothing 
length before the actual shock front (which travels to the right in 
Fig. II It. whereas the M&M switch turns on viscosity later, lagging 
our method by about four particle separations. As a consequence, 
the transition of the fluid values across the shock front is slightly 
smoother with our method than with the M&M method. 

Note that the irregularities around the contact discontinuity at 
J = 0. 138 common to all schemes tested are not related to artificial 
viscosity (the irregularities in a at that point could be removed by 
choosing non-zero initial a at the in itial discontinuity) ; they can be 
alleviated by artificial conductivity ( |Pricel2004l200i) . 
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Figure 13. Particle positions in tlie x-y plane of 3D simulations of a A1 = 
20 ram shock along x direction. Particles are coloured red if there initial 
positions was ,ro < -0.45 and green if xo > -0.45 



Figure 14. Shearing shock test: density and velocity for various SPH 
schemes (symbols) and a grid-code simulation (curve). Initial velocities are 
given by equation )2U with s = lOOf = 100c. 



5.2 Strong shocks and particle penetration 

In 33.4l and Fig. [3] we already demonstrated that our new method is 
superior to the M&M scheme in resolving 'subsonic shocks' (ve- 
locity discontinuities smaller than the sound speed) and compara- 
ble in resolving shocks of Mach number ~ 1. Here, we extend this 
comparison to high Mach numbers. Fig. ll2| shows the result for the 
ID ram test with M = 50. The lPricei ilOOj) version of the M&M 
method, which uses equation l IlO) with ormax = 2, is implemented in 
some contemporary SPH codes, and has been used in our tests so 
far, fails this test: a remains too low and as a consequence the ve- 
locity discontinuity is not correctl y smoothed and some post-s hock 
ringing occurs. To give credit to iMorris & Monagliaril ( Il997h , we 
also tested their original method and find it to work well (stars in 
Fig.inj. Our new method works about as well as the original M&M 
scheme, with a reaching the same level, though our scheme detects 
the coming shock much earlier: a is ahead of the original M&M 
method by about four particle separations. 

Whilst the main role of artificial viscosity is to resolve shocks 
by transferring entropy, a secondary but vital role is to prevent 
particle penetra tion, which requires strong viscosity in high Mach 
number shocks . Ib atd l l 1 9951 performed many tests to determine the 
value of the parameters a and /3 needed to prevent particle pene- 
tration in ram shock tests of various Mach numbers. For particles 
arranged in face-centred-cubic or cubic grids. Bate found that ap- 
propriate values for the viscosity parameters can prevent particle 
penetration for shock s up to = 8. Most SPH practitioners opt 
for a value of /? = 2a |m orris &Monaghanll 19971) . 

To determine the correct value of fS required for the new 
scheme, we perform high-resolution 3D runs of ram shocks with 
M = 20 and various values for 0/a. We smooth th e initial veloc- 
ity discontinuity, as suggested bv lMonaghai] h997l) . to provide the 
method with a situation realistic for SPH, such as would have arisen 
for a shock forming from continuous initial conditions. 

For different values of /3/a with our viscosity scheme and the 
two variants of the M&M switch, we plot in Fig. [T3] the x and y 
positions (for all values of z) of particles near the shock front at a 
late time. The colour coding distinguishes particles which at that 
time should be up- (red) or downstream (green). Our scheme pre- 



vents particle penetration with /3 = or (for /3 = there is one layer 
of overlap). The original M&M scheme with the st andard choice 
/3 = 2a also avoids particle penetration, but not the I Pried i2004) 
version, again a consequence of too little viscosity. 



5.3 A shearing shock 

This test combines a shock with a perpendicular shear and presents 
a difficult test for any SPH scheme. We use periodic boundary con- 
ditions and start from a face-centred cubic grid and velocities 

Ui = -fosign(x), Vy = s sm(nx), and v- = 0. (21) 

In Fig. [141 we present results for various SPH simulations as well 
as a grid-code run for s = I00i5li = 1 00c. The M&M method pro- 
duces a large viscosity due to the shear-induced errors in V-v, lead- 
ing to spurious results. Using the Balsara limiter with either M&M 
or Standard SPH gives in much better results, though the shock is 
clearly over-smoothed. The new scheme is able to limit the viscos- 
ity to the correct level, allowing good capturing of the shock and 
retaining particle order in the post-shock region. 

Note that this is a difficult test for any SPH implementation: 
without viscosity reduction (as in standard SPH) the shear flow is 
strongly damped, while viscosity reduction schemes (M&M as well 
as ours) suffer from the problem of shear-induced errors. These po- 
tentially result in too much viscosity and over-smoothing of the 
shock. Our limiter was able to control this problem, but for yet 
larger ratios s/6v of shear to shock amplitude this problem becomes 
too difficult for any SPH implementation. 



5.4 Evrard Test 

In this test the inward gravitational pull of a gas cloud exceeds its 
outward pressure force causing the cloud to collapse under its own 
self-gravity. The initial conditions consist of a gas sphere with den- 
sity profile (iEvrard.1988.) 

M I 
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Figure 15. The Evrard test (see text for the initial setup): shown are various physical quantities (K = Pp^'' is the entropy function) and a at different times for 
SPH simulations wit h N = particles using our new viscosity scheme (blue) or the original M&M method (red). Also shown (black) are the results from 
ID PPM calculation jSteinmetz & Mullej|l993h . Not every particle is plotted. 



for r <R and p = Q for r > R. Initially the gas is at rest and has con- 
stant specific internal energy u = 0.05 GM/R, which corresponds to 
a virial ratio -2U/W = 0.075 ^ 1. The initial gravitational inward 
pull is the same at each radius, while the pressure forces decline 
outwards, leading to collapse and, as a consequence, formation of a 
shock, which steepens and evolves into a strong shock propagating 
outwards as more incoming material joins the jam. Even though the 
problem is initially spherically symmetric, the SPH realisation of 
initial conditions cannot be exactly spherically symmetric and the 
system may well evolve away from sphericity, for instance driven 
by dynamical instabilities. 

We use a unit system such that G = R = M = 1 and represent the 
cloud by 100280 SPH particles, initially placed on a face-centred- 
cubic grid which is then radially stretched to match the density. 
Fig.[T5]compares the simulation results for our method, the origina l 
M&M method, and a ID calculation bv lSteinmetz & Miillerl n993h 
using the piece-wise parabolic method (PPM). 

At early times (t = 0.39, left column) the results from all three 
methods match very well, but the M&M scheme already shows 
a large viscosity. At later times a shock forms (at r x 0.13 by 
t = 0.78), which moves outwards until it reaches the end of the 
sphere, when a significant fraction of the gas still has outwards 



velocities (by t = 1.95). The most obvious difference between 
the two SPH schemes is the amount of (artificial) dissipation: the 
M&M method is much more viscous, resulting in significant over- 
smoothing of the shock front by ? = 0.78 accompanied by unphysi- 
cal pre-shock heating as visible in the entropy (K) profile. Our new 
scheme agrees better with the ID calculation, in particular in the 
inner (post-shock) regions. Note that with our new method a peaks 
well before the shock arrives (at / = 1.17), while for the M&M 
method the peak in a appears actually slightly after the shock. 

We found this a valuable test as early versions of our scheme 
tended to be far too viscous, while our final version passes this test 
ahead of the M&M switch. Standard SPH (not shown in the figure) 
shows similar results, though the shock at ? = 0.78 appears less 
smoothed than with the M&M method but more smoothed than 
with the new scheme. 



6 SUMMARY 

Any hydrodynamical numerical method requires some form of ar- 
tificial viscosity in order to resolve shocks (in grid methods, ar tifi- 
cial viscosity is implicit in the Riemann solver. lMonaghanI 19971) . In 
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grid codes, such as Ramses ( lTevssiej200"3) . interpolation methods 
are employed to effectively suppress artificial viscosity away from 
shocks. Most SPH simulations to this date hardly use such precau- 
tions and, as a consequence, adiabatic oscillations and shear-flows 
are damped. Note that this affects state-of-the-a rt simula ti ons of , 
e.g. galaxy formation, which usually only emplov lSalsaral 's ( Il995h 
rather inefficient method to reduce some adverse effects of artificial 
viscosity on rotation discs^ 

The method of Morris & MonaghanI jl997h . which reduces the 
default amount of artificial viscosity by an order of magnitude com- 
pared to standard SPH practice, has only recently been recognised 
as advantageous. In this method, explained in detail in 92.31 in- 
dividual artificial viscosities a, are adapted by integrating a dif- 
ferential equation. Though constituting a major improvement, this 
method remains unsatisfactory, because it still damps adiabatic os- 
cillations and over-smoothes weak shocks, as we argued in ^and 
demonstrated in ^ 

In we present a novel method, which improves upon that 
of|M orris & MonaghanI in four important ways. 

• We set ffn,i„ = enabling cy, — > away from shocks and effec- 
tively modelling the fluid as inviscid. 

• We use V-y = d(V-y)/d/ < rather than V-u < as shock 
indicator. This distinguishes pre-shock from post-shock regions 
(where V-u > but V u < 0) and discriminates much better be- 
tween converging flows and weak shocks. 

• We set or, directly to an appropriate local value a^^c, instead of 
growing it by integrating a differential equation. 

• We use an improved estimator for V-u and V-u and employ a 
limiter to avoid viscosity driven by shear-induced errors. 

Together these novelties result in a significantly improved artificial 
viscosity method, in particular the viscosity is increased to an ap- 
propriate level well before an incoming shock. The implementation 
details, i.e. the precise way of setting aioc from V-u and the ex- 
act form of the limiter, may well be subject to improvements. Any 
reader who considers modifying these details is advised to consider 
the behaviour of the resulting method for a test suite comprising 
noise suppression as well as shear and strong shocks, for example 
the tests of Figures g] (8] and[T4l 

For static equilibria V w = and i> = 0, and our new shock in- 
dicator (as well as the M&M shock indicator) are only triggered by 
velocity noise. As long as particle order is maintained, such noise 
triggers only negligible amounts of viscosity, unlike the situation 
with the M&M method, whose minimum viscosity a^^^^ = O.I is 
often sufficient to affect the simulations (as demonstrated in ®. 
Nonetheless, the noise-induced viscosity is sufficient to suppress 
particle disorder, as demonstrated in 33.51 

For dynamic equilibria V u = (and V-u = 0) but v t 0. 
However, in multi-dimensional flows strong shear induces false de- 
tections of V-i/ (and V-w), even with best possible particle order, 
for reasons explained in Appendix IB II In simulations of differen- 
tially rotating discs, this problem strongly affects the M&M method 
(even when using the Balsara switch). We avoid this problem by 
applying a limiter (see 33.3! ) as well as using improved estimators 
for V v and V-w, see Appendix IB2l for details. (Alternatively, if no 
strong shear flows are present, the standard estimators should suf- 
fice, though still in conjunction with a limiter using |V x u| as a 
proxy for the shear amplitude.) 

These improved estimators also provide the full velocity and 
acceleration gradient matrices for each particle (and increase the 
computational costs by ~ 30%). The knowledge of the velocity gra- 



dient matrix V and its traceless symmetric part, the shear S, is also 
useful for implementing physical viscosity 

pi/ = V-[77S + ftr(V)] (23) 

(with rj and f the shear and bulk viscosity coefficients) in SPH. 

In sections [331 l4l and|5] we demonstrate convincingly that our 
technique successfully deals with the following four situations. 

Shocks are resolved at least as well, if not better, than with any 
previous technique; 

adiabatic oscillations, such as sound waves or stellar pulsations, 
remain undamped over many periods, which was not possible with 
any previous SPH implementation; 

strong shear flows, such as in accretion discs, are modelled vir- 
tually inviscid, while shearing shocks are well resolved without be- 
ing over-smoothed; 

particle disorder is suppressed at least as well as with the M&M 
method. 

In particular, in the regime of convergent flows and weak shocks 
our new method is far superior to any previous scheme, which all 
required a significant increase in resolution just to suppress adverse 
effects of artificial viscosity. 
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APPENDIX A: DETAILS OF THE SPH SCHEME 

For completeness, we give here a brief description of our SPH 
method, which is largely similar to previous methods, but may dif- 
fer in some details. 



Al Density and adaptive smootliing lengths 

Let V denote the number of spatial dimensions, then we adapt the 
individual smoothing lengths /j, such that p, = M;, with M;, = 
inNh/Vv a global constant, defined in terms of the number Ni, of 
neighbours, the mass m of each SPH particle, and the volume V,, of 
the unit sphere. In this work, we use N/, = 5, 13, and 40 for v = 
1, 2, and 3 dimensions, respectively. Inserting the density estimator 
(TJ, we find 



h-Pi = ^ mj w(nj) with r,y = 



(Al) 



where we have re-written the SPH kernel as Wdx^l, /?,) = hj^w(rij) 
with the dimensionless fu nction w(r). For this work, we employ the 
usual cubic spline kernel jMonaghan & Lattanzid[l983) 



w(r) ■■ 




■r) r< 1/2, 

1/2 < r < 1, (A2) 
otherwise. 



At each time step, the A, are adjusted by performing one Newton- 
Raphson step in log h-\og(h^p) space, i.e. 



hi «- hi 



Ml, 



fi/v 



with a factor of order unity 

Y^jmjWij 



fi = -v 



ZjtnjrjjWij 



(A3) 



(A4) 



where w,j = w(rij) and w(r) = w'(r)/r. This method converges 
extremely well, except when /i, was much too small. In this 
case, faster convergence can be achieved by subtracting the self- 
contribution (which does not depend on hi). Thus, whenever /z!p, < 
Ml, we use instead of l lA3b 

^_ ^ ^_ /M,.-m,w(0)V^> ^. ^ - ^ ^^3^ 
\ h]pi - nii w(0) / 2 m, rf. Wij 



The time derivatives hi are obtained by demanding d(/jj'p,)/df = 0: 
hi 2 j nij Vij ■ Xij Wij 



hi Y^jinjxlwij 



A2 Pressure forces 



(A6) 



We use SPH equations of motion derived from the simple SPH 
Lagrangian £ = Y.kmkikxl - ut 
di(/dp = P/p', this gives 



). Together with the relatioiQ 



Vi=- 



nii dXi 



P.f. 



"21 v+2 ^ -7, F+2 



(Al) 



where the factors / and fj (eauation lA4t arise from the fact that the 
derivatives dpkldxi have to be taken at fixed /i^p^.. The work done 
by these pressure forces has to be balanced by 



P,k P.f, 



^nijUij - Xij Wij. 



P. hi pjhy^ ^ 



A3 Artificial viscosity 

For the artificial viscosity drag and heating we actually use 



(A8) 



H,^ 

(y,)AV = - ^jOTyXy — 

J 

(",)aV = ^ inj Vij ■ X,; 



aafi - 



P.h] 
2 p,/,-^^" 



(A9) 
(AlO) 



with Hy = -pij(cij - b^ij), where 



b 



for Vij ■ Xij < 0, 



otherwise; 



(All) 



while the parameter b has the meaning of pi a for traditional SPH. 
Note that equations l |A9t and l [2at differ only by 0(hfj). The differ- 
ence between equations dAlOt and l l2bb is more pronounced since, 
similarly to equation jASi , we do not symmetrise the contributions 
w.r.t. / and j. 



A4 Time Integration 

Our scheme employs a kick-drift-kick leap-frog time integrator, 
which is second-order accurate. With this scheme, a full (global) 
time step of size 6t consists of the following sub-steps ('«— ' means 
'is replaced by'). 



initial kick Compute i», and m, at half step 

£), = Vi + J dtVi, 
Ui = Ui + jStiii. 

full drift Advance / and Xi by a full step: 

t (- t + 6t 
Xi «— Xi + Stiii. 



(A12) 



(A13) 



* Alternatively, for an ideal-gas equation of state one may replace u in the 
Lagrangian with u = Kp^^^ /{y - 1) and consider the entropy function K = 
Pp^''' to be constant iSpringel & Hemquislll20o3) . 
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prediction Predict W;, m,, and A, at full step: 



- w, + 6tVi, 

-;(, exp((5f «,/«,), (A14) 
hj <— hjexp(6thj/hj). 

sweep Compute h]'pj and fi (equations lAll and lA4b . 
adapt Adjust (equation |A3| or |A5t . 

sweep 1 Compute p,, and hj (eqs. lAlllA4l and lA6t as well as 
V w,, V vj, and Rj (eqs. lB8IIB12l andll7l using v and V v from the 
previous time step). 

between sweeps Obtain P, and c, from p, and w, via the equation 
of state, and adapt a, via (using eQS.[9landll4b 



a, 



J ffloc if ff, < Q-ioc 

1 ttioc + (o'; - Cloc) exp(-5f/T,) Otherwise. 



(A15) 



sweep 2 Compute i/, (eqs. IA7l and |A9| plus gravitational forces) 
and iij (eqs. |A8| and |A10| plus external heating or cooling), 
final liiclt Set w, and at full step: 



v, = Vj + J 6tVj, 
iij = iij + ^ Sliij. 



(A 16) 



In the initial kick and prediction steps, the time derivatives are 
known from the previous time step (in case of the very first time 
step, they need to be precomputed). Note that the quantities pre- 
dicted in l lA14t enter the final u, and m, only indirectly via the com- 
putation of the time derivatives. 

We use an oct-tree, generated just before sweep 0, to find all 
interacting particle pairs, which are then remembered in an inter- 
action list, whereby allowing for the fact that hj may grow slightly 
during adjustment (just after sweep 0). Utilising this interaction list 
in sweeps 1 and 2 is much faster than further tree walks. The same 
oc t-tree is also us ed in computing gravitational forces, as outlined 
bv lDehnenlJiooih . 

Our scheme can also be implemented with adaptive individual 
time steps organised in a hierarchical block-step scheme, though 
we have not used this in the tests presented in this study. 



APPENDIX B: ESTIMATING V y AND V v 

Bl Failure of the standard SPH estimator for V w 



Our constraint that /ij'p, be constant (see gAT) implies p,7p, = 
-vhj/hj. Together with the continuity equation p + pV-v = and 
equation l lA6t this yields the simple velocity-divergence estimate 



V-w, = V- 



(Bl) 



While this estimate satisfies the continuity equation for the SPH 
density estimate p,-, it is not necessarily accurate. To see this, con- 
sider the matrix (0 denotes the outer or dyadic vector product) 



(B2) 



with Wjj some weighting factor. Assuming a smooth velocity field, 
we may replace Vy in equation l lB2t with its Taylor expansion Vy = 
V, • Xy + 0(\xy\^), where V, = V (8> is the gradient of v at position 
Xi, and obtain 



D, = V,-T + h.o.t. 

with the symmetric matrix 

Tj = 'Zj Xy Xy Wy. 



(B3) 
(B4) 



Comparing l lB2t and l lB4t to the simple estimator JBU . we see 
that the latter corresponds to (conveniently dropping the index i) 
V v = vtr(D)/tr(T) and the weights ivy = mjWy. If we split V into its 
isotropic part (divergence), the symmetric traceless part S (shear), 
and the antisymmetric part R (vorticity). 



Vyl + S + R, 



(B5) 



and insert it into iB3i , we find for the simple estimator dBU 

= V-u + vtr(S-f)/tr(T) + h.o.t. (B6) 

where T denotes the anisotropic (traceless) part of T. Thus, the sim- 
ple estimator JBU contains an 0(h°) error term, which originates 
from anisotropy of T in conjunction with velocity shear (owing to 
the symmetry of T the vorticity is harmless). For perfectly symmet- 
ric particle distributions T = 0, but in general T # such that in the 
presence of strong shear even a small residual T results in a fail- 
ure of the simple estimator JBU . This typically happens in differ- 
entially rotating discs, where (i) the velocity field is divergent-free 
but contains shear and (ii) even in the absence of noise J t owing 
to the shearing particle distribution. 

B2 A more accurate V-y estimator 



From equation iB3i , we can also estimate 
V, = D, I" ', 

which allows an improved divergence estimator 
V^, = tr(D,-Tr'). 



(B7) 



(B8) 



In order to assess the error of this estimator, let us expand the flow 
to second order, replacing equation iB3\ with (dropping the index ; 
and using suffix instead of matrix notation) 



with the symmetric tensor U, 
into l lB7t we find 



(B9) 

YijXy Xy Xy Wy. Inserting this 



UvA„T„„ -I- h.o.t.. 



(BIO) 



Thus, while this estimator avoids an Oih*^) error, we still have an 
(?(/j') error term (since U is one order higher in h than T). We can 
reduce the 0{h^) error by a careful choice of the weights wy. If, 
for instance, wy = mjwy/pj then U — » to leading order in the 
continuum limit by virtue of the isotropy of the kernel. This limit, 
which is commonly used to assess SPH estimators, replaces Yij t^j 
with f p(Xj) dXj under the assumption of a smooth density without 
particle nois^H As these conditions are hardly ever truly satisfied, 
we can only reduce but not eliminate the 0(h^) error term — as we 
do not even try to avoid the 0(h^) error (hidden in 'h.o.t.' above), 
such a reduction should be okay in most cases. 



B3 Estimating V-u 

We can estimate V-y either from the change in the estimated V-w 
over the last time step or as the trace of V, the total time derivative 
of V. Since (with A = V w the gradient of the acceleration) 



V = A - V 



(BID 



^ Under these conditions also T, which causes the O(h^) error term in the 
simple V-K estimator, vanishes. 
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(a good exercise for your undergraduate students), we can estimate 

= tr(A, - V"). (B12) 

Here, the estimate A, is obtained from the accelerations at the pre- 
vious time step in a way analogous to estimating V,, in particular 
we need to compute the matrix X and its inverse only once. The 
lowest-order error in this estimate again is 0(h^) oc U,, such that 
reducing U, by careful choice of the weights remains a good idea. 

Note that, by virtue of equation dBllb . we could estimate V-i> 
also as V-i> - tr(V") with the acceleration divergence V-u estimated 
using the standard divergence estimator, in the hope that its 0{h°) 
error term is small since the acceleration is hardly sheared. 



